using LinearAlgebra, Plots, Printf, SparseArrays, LaTeXStrings
We first discuss how to solve
\begin{align*} u_t + \alpha u_{x} &= \kappa u_{xx}, \quad a < x < b, \quad t > 0,\\ u(a,t) &= g_0(t),\\ u_x(b,t) &= 0,\\ u(x,0) &= \eta(x). \end{align*}One approach would be to combine everything and use one (implicit) method to time step the MOL discretization. Here one might want to use something like TR-BDF-2 (as opposed to trapezoid) so that if the eigenvalues from advection that lie on the imaginary axis leak into the right-half plane, after being perturbed, the method is still most likely stable. Although trapezoid is probably sufficient.
To illustrate splitting methods, we take a different approach and treat each term separately. In this lecture we won't do Strang splitting. We will just use the first-order (in time) approach (11.17).
a = 0.; b = 10
h = 0.01
m = convert(Int64,b/h) #increase by 1 over Dirichlet conditions
T = 4
κ = .1;
k = h
B = Tridiagonal(fill(1.0,m-1),fill(-2.0,m),fill(1.0,m-1))
B[end,end-1] = 2
B *= 1/(h^2)
B1 = I - κ*(k/2)*B
B2 = I + κ*(k/2)*B;
r = κ*k/(2*h^2)
g0 = t -> 2*sin(2*pi*t).^2
η = x -> 0.
#3 (generic function with 1 method)
n = convert(Int64,ceil(T/k))
x = h:h:b |> Array
U = map(η,x)
t = 0.0
anim = Animation()
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
# diffuse with Crank-Nicolson
U = (B2*U)
U[1] += r*(g0(t)+g0(t-k))
U = B1\U
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
end
end
gif(anim,"no_advection.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/no_advection.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
Recall the setting of Lax-Friedrichs:
$$U^{n+1} = U^n - k A U^n + \frac{h^2}{2} B U^n$$where $B$ is the second-order diffusion matrix. But we already have a matrix set up for diffusion in our problem so we could use that.
There is a problem though. As the advection parameter $\alpha$ increases, our required step size $k \leq h/|a|$ decreases relative to $h$. And as $k$ decreases relative to $h$, the effect of the diffusion term $h^2/2B$ will become more and more pronounced. We then have to decrease $h$ further to maintain the same accuracy. We have an eye toward going to two spatial dimensions and this will be way too costly. So, we need to use the Lax-Wendroff approach.
Recall the derivation of the Lax-Wendroff method. If our MOL system with boundary conditions is given by
\begin{align*} U'(t) = A U(t) + g(t). \end{align*}Here $A = - \frac{\alpha}{2h}S$ where $S$ is a matrix with $-1$'s on the first subdiagonal and $1$'s on the first superdiagonal. And $g(t) = \frac{\alpha}{2h}g_0(t) e_1$
To derive Lax-Wendroff with a Dirichlet condition at the left, we first ignore $g$
\begin{align*} U(t+k) &= U(t) + k U'(t) + \frac{k^2}{2} U''(t) + O(k^3)\\ &= U(t) + k A U(t) + \frac{k^2}{2} A U'(t) + O(k^3)\\ & = U(t) + k A U(t) + \frac{k^2\alpha^2}{8h^2} S^2 U(t) + O(k^3). \end{align*}Then the key step in the Lax-Wendroff derivation is to replace $\frac{1}{4h^2}S^2$ with $B$, giving
\begin{align*} U^{n+1} = U^n + k A U^n + \frac{k^2 \alpha^2}{2} B U^n. \end{align*}At a typical grid point this becomes
\begin{align*} U_j^{n+1} = U_j^n - \frac{\alpha k}{2 h} \left( U_{j+1}^n - U_{j-1}^n \right) + \frac{\alpha^2 k^2}{2 h^2} \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n \right). \end{align*}Evaluating this at $j =1$ and using $U_0^{n} = g_0(t_n)$ tells us how to impose the boundary conditions.
a = 0.; b = 10
h = 0.01/2
m = convert(Int64,b/h) #increase by 1 over Dirichlet conditions
T = 4
α = 3.
κ = .1;
k = h/(2α)
B = Tridiagonal(fill(1.0,m-1),fill(-2.0,m),fill(1.0,m-1))
B[end,end-1] = 2
B *= 1/(h^2)
B1 = I - κ*(k/2)*B
B2 = I + κ*(k/2)*B;
r = κ*k/(2*h^2)
A = Tridiagonal(fill(-1.0,m-1),fill(0.0,m),fill(1.0,m-1))
A[end,end-1] = 0
A *= -α/(2h);
A1 = I + k*A + (k^2*α^2)/(2)*B; # Lax-Wendroff matrix
g0 = t -> 2*sin(2*pi*t).^2
η = x -> 0.
#7 (generic function with 1 method)
n = convert(Int64,ceil(T/k))
x = h:h:b |> Array
U = map(η,x)
t = 0.0
anim = Animation()
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
## O(k) splitting method
# diffuse with Crank-Nicolson
U = (B2*U)
U[1] += r*(g0(t)+g0(t-k))
U = B1\U
#advect with Lax-Wendroff
U = A1*U
U[1] += (α*k/(2h) + (α^2*k^2/(2h^2)))*g0(t-k)
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
end
end
gif(anim,"advection_diffusion.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_diffusion.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
We then can simulate a rudementary "wave tank" by decreasing the diffusion and turning off the boundary condition after some time.
a = 0.; b = 10
h = 0.01
m = convert(Int64,b/h) #increase by 1 over Dirichlet conditions
T = 6
α = 3.
κ = 0.01;
k = h/(2α)
B = Tridiagonal(fill(1.0,m-1),fill(-2.0,m),fill(1.0,m-1))
B[end,end-1] = 2
B *= 1/(h^2)
B1 = I - κ*(k/2)*B
B2 = I + κ*(k/2)*B;
r = κ*k/(2*h^2)
A = Tridiagonal(fill(-1.0,m-1),fill(0.0,m),fill(1.0,m-1))
A[end,end-1] = 0
A *= -α/(2h);
A1 = I + k*A + α^2*k^2/(2)*B; # Lax-Wendroff matrix
g0 = t -> t < .5 ? 2*sin(2*pi*t).^2 : 0.0
η = x -> 0.
#11 (generic function with 1 method)
n = convert(Int64,ceil(T/k))
x = h:h:b |> Array
U = map(η,x)
t = 0.0
anim = Animation()
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
## O(k) splitting method
# diffuse with Crank-Nicolson
U = (B2*U)
U[1] += r*(g0(t)+g0(t-k))
U = B1\U
#advect with Lax-Wendroff
U = A1*U
U[1] += (α*k/(2h) + (α^2*k^2/(2h^2)))*g0(t-k)
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
end
end
gif(anim,"linear_wavetank.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/linear_wavetank.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
We now discuss how to solve the nonlinear problem
\begin{align*} u_t + \alpha u u_{x} &= \kappa u_{xx}, \quad a < x < b, \quad t > 0,\\ u(a,t) &= g_0(t),\\ u_x(b,t) &= 0,\\ u(x,0) &= \eta(x). \end{align*}Now, if we tried to combine everything into a single implicit method, we would have to solve a nonlinear system of equations at each time step. The splitting approach avoids this. Handling the diffusion is therefore the same. We have to think more carefully about how to handle the nonlinear advection term. We want a Lax-Wendroff-like approach but just replacing $\alpha$ with $\alpha u$ does not give a second-order method.
We follow the approach just below (10.18) in the text and just Taylor expand:
\begin{align*} u_t &= - \alpha u u_x,\\ u_{tt} & = - \alpha u_t u_x - \alpha u u_{xt}\\ & = 2\alpha^2 u u_x^2 + \alpha^2 u^2 u_{xx}\\ u(x,t+k) &= u(x,t) + k u_t (x,t) + \frac{k^2}{2} u_{tt}(x,t) + O(k^3),\\ u(x,t+k) &= u(x,t) - \alpha k u(x,t) u_x (x,t) + \frac{k^2}{2} \left( 2\alpha^2 u(x,t) (u_x(x,t))^2 + \alpha^2 (u(x,t))^2 u_{xx}(x,t)\right)+ O(k^3) \end{align*}We now replace derivatives with centered differences:
\begin{align*} U^{n+1}_j = U^n_j - \frac{\alpha k}{2h} U_j^n \left( U_{j+1}^n - U_{j-1}^n \right) + \frac{\alpha^2 k^2}{4h^2} U_j^n \left( U_{j+1}^n - U_{j-1}^n \right)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_j)^2 \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n \right). \end{align*}Setting $j =1$ gives the relation for the boundary:
\begin{align*} U^{n+1}_1 &= U^n_1 - \frac{\alpha k}{2h} U_1^n \left( U_{2}^n - g_0(t_n) \right) + \frac{\alpha^2 k^2}{4h^2} U_1^n \left(U_{2}^n - g_0(t_n) \right)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 \left( g_0(t_n) - 2 U_1^n + U_{2}^n \right) \end{align*}We write this as \begin{align*} U^{n+1}_1 &= F(U_j^n,U_{j+1}^n) + \frac{\alpha k}{2h} U_1^ng_0(t_n) - \frac{\alpha^2 k^2}{2h^2} U_1^n U_2^n g_0(t_n) + \frac{\alpha^2 k^2}{4 h^2}U_1^ng_0(t_n)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 g_0(t_n). \end{align*} or in matrix-vector form: \begin{align*} U^{n+1} &= U^n - k\,\mathrm{diag}(U^n) (A U^n) + k^2 \mathrm{diag}(U^n) (A U^n)^2 + \frac{\alpha^2 k^2}{2}\mathrm{diag}(U^n)^2 (B U^n) + h^n e_1,\\ h^n &= \frac{\alpha k}{2h} U_1^ng_0(t_n) - \frac{\alpha^2 k^2}{2h^2} U_1^n U_2^n g_0(t_n) + \frac{\alpha^2 k^2}{4 h^2}U_1^ng_0(t_n)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 g_0(t_n). \end{align*}
a = 0.; b = 5
h = 0.01
m = convert(Int64,b/h) #increase by 1 over Dirichlet conditions
T = 6
α = 2.
κ = 0.01;
k = h/4
B = Tridiagonal(fill(1.0,m-1),fill(-2.0,m),fill(1.0,m-1))
B[end,end-1] = 2
B *= 1/(h^2)
B1 = I - κ*(k/2)*B
B2 = I + κ*(k/2)*B;
r = κ*k/(2*h^2)
A = Tridiagonal(fill(-1.0,m-1),fill(0.0,m),fill(1.0,m-1))
A[end,end-1] = 0
A *= -α/(2h);
g0 = t -> t < 1. ? 2*sin(2*pi*t)^2 : 0.0
η = x -> 0.
#15 (generic function with 1 method)
n = convert(Int64,ceil(T/k))
x = h:h:b |> Array
U = map(η,x)
t = 0.0
anim = Animation()
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
# diffusion
U = (B2*U)
U[1] += r*(g0(t)+g0(t-k))
U = B1\U
# nonlinear advection
AU = A*U
U1 = U[1]
U2 = U[2]
U = U + k*U.*AU + k^2*U.*(AU.^2) + (α^2*k^2/2)*(U.^2).*(B*U)
g = g0(t-k)
H = (α*k/(2h))*U1*g - (α^2*k^2/(2h^2))*U1*U2*g + (α^2*k^2/(4h^2))*U1*g^2 + (α^2*k^2/(2h^2))*U1^2*g
U[1] += H
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [a,b], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([a,b],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
end
end
gif(anim,"nonlinear_wavetank.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/nonlinear_wavetank.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
ay = 0.; by = 1
hy = 0.01
my = convert(Int64,(by-ay)/hy) + 1 #increase by 2 over Dirichlet conditions
T = 2
α = 2.
κ = 0.01;
k = hy/4
By = Tridiagonal(fill(1.0,my-1),fill(-2.0,my),fill(1.0,my-1))
By[end,end-1] = 2
By[1,2] = 2
By *= 1/(hy^2)
B1y = I - κ*(k/2)*By
B2y = I + κ*(k/2)*By;
ry = κ*k/(2*hy^2)
Ay = Tridiagonal(fill(-1.0,my-1),fill(0.0,my),fill(1.0,my-1))
Ay[1,2] = 0
Ay[end,end-1] = 0
Ay *= -α/(2hy);
η = x -> exp(-25(x-.5)^2)
#17 (generic function with 1 method)
n = convert(Int64,ceil(T/k))
y = ay:hy:by |> Array
U = map(η,y)
t = 0.0
anim = Animation()
plot(y, U, xaxis = [ay,by], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([ay,by],[U[1],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B1y\(B2y*U)
AU = Ay*U
U = U + k*U.*AU + k^2*U.*(AU.^2) + (α^2*k^2/2)*(U.^2).*(By*U)
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(y, U, xaxis = [ay,by], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t), fill = (-2,:lightblue))
plot!([ay,by],[U[1],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
frame(anim)
end
end
gif(anim,"burgers_neumann.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/burgers_neumann.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
ax = 0.0; bx = 5.0;
ay = 0.0; by = 1.0;
hx = 0.1/10
hy = 0.1/2
x = ax+hx:hx:bx |> Array
y = ay:hy:by |> Array
X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(reverse(y), 1, length(x));
Y = Y[end:-1:1,:]
αx = 2.0
αy = .1
κ = 0.01
k = min(hx,hy)/(4*max(αx,αy))
my = convert(Int64,(by-ay)/hy) + 1 #increase by 2 over Dirichlet conditions
By = Tridiagonal(fill(1.0,my-1),fill(-2.0,my),fill(1.0,my-1))
By[end,end-1] = 2
By[1,2] = 2
By *= 1/(hy^2)
B1y = I - κ*(k/2)*By
B2y = I + κ*(k/2)*By;
ry = κ*k/(2*hy^2)
Ay = Tridiagonal(fill(-1.0,my-1),fill(0.0,my),fill(1.0,my-1))
Ay[1,2] = 0
Ay[end,end-1] = 0
Ay *= -αy/(2hy);
mx = convert(Int64,(bx-ax)/hx) #increase by 1 over Dirichlet conditions
Bx = Tridiagonal(fill(1.0,mx-1),fill(-2.0,mx),fill(1.0,mx-1))
Bx[end,end-1] = 2
Bx *= 1/(hx^2)
B1x = I - κ*(k/2)*Bx
B2x = I + κ*(k/2)*Bx;
rx = κ*k/(2*hx^2)
Ax = Tridiagonal(fill(-1.0,mx-1),fill(0.0,mx),fill(1.0,mx-1))
Ax[end,end-1] = 0
Ax *= -αx/(2hx);
function plot_wave(U,x,y,t,cl,width=800)
l = @layout [a;b]
p1 = surface(x, y, U, camera = (10,70), zaxis = [cl[1],cl[2]], clims= cl, aspectratio = 1.5, xlabel = L"x", ylabel = L"y", zlabel = L"u(x,y,t)")
p2 = contour(x, y, U, clims=cl, fill = true, aspectratio = 1, xlabel = L"x", ylabel = L"y")
plot(p1, p2, layout = (2,1), size = (width, 7*width/10), title = @sprintf("t = %1.4f",t))
end
plot_wave (generic function with 2 methods)
η = (x,y) -> 0*exp(-25*(x-2)^2-50*(y-.5)^2)
g0 = (y,t) -> 2*sin(2*pi*t).^2 .*exp(-25*(y-.5)^2)
U = map(η,X,Y)
T = 10
n = convert(Int64,ceil(T/k))
t = 0.0
cl = (-1,2)
plot_wave(U,x,y,t,cl) |> IJulia.display
anim = Animation()
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
# diffuse in y
U = B1y\(B2y*U)
# advect in y
AU = Ay*U
U = U + k*U.*AU + k^2*U.*(AU.^2) + (αy^2*k^2/2)*(U.^2).*(By*U)
U = U' |> Array
# diffuse in x
U = (B2x*U)
U[1,:] += map(yy -> rx*(g0(yy,t)+g0(yy,t-k)),y)
U = B1x\U
# advect in x
AU = Ax*U
U1 = U[1,:]
U2 = U[2,:]
U = U + k*U.*AU + k^2*U.*(AU.^2) + (αx^2*k^2/2)*(U.^2).*(Bx*U)
g = map(yy -> g0(yy,t-k), y)
H = (αx*k/(2h))*U1.*g - (αx^2*k^2/(2h^2))*U1.*U2.*g + (αx^2*k^2/(4h^2))*U1.*g.^2 + (αx^2*k^2/(2h^2))*U1.^2 .*g
U[1,:] += H
U = U' |> Array
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot_wave(U,x,y,t,cl) |> IJulia.display
frame(anim)
end
end
gif(anim,"viscous_Burgers_2D.gif")